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I, ABSTRACT 


A method of calculating the eigenvalues and the expansion coefficients 
of sphero:dal wave functicns for the prolate and oblate case using a high- 
speed computer is presented. Special emphasis is placed on a new explanation 
of the role the parameter n plays mm the numerical determination of eigen- 
values. Curves and tables of eigenvalues are included. 
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il. INTRODUCTION 


Spheroidal wave functious are important im numerous fielas of physics, 
especially in the field of acoustics. A greater interest in these functions 
has brought about an ancreasea demand for more elaborate tabuiation of 
values of these functions and quantities related to them. 

In answer to this demand the Theory Division of the Applied Mathematics 
Laboratory has completed a study leadmg to the numerical determimation of 
spheroidal eigenvalues and expansion Coefficients of the prolate and oblate 
case. In our nomenclature we follow the notation of Flammer?, but include 
work in ranges of m c, andn exceeding those of Flammer. 

Work in these extended ranges was made possible by a re-examination 
of the parameter n. The parameter n 1s discussed at length in the following 
sections. 
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Ill) NUMERICAL DETERMINATION OF EIGENVALUES 


Eigenvalue Equations 
In his book ''Spheroidal Wave Functions, ' Flammer (2, p. 16 | studies 


the differential equation satisfied by the angle furctions Sant’? n), where 


in his notation 


dS. -c, ni 2 
d ‘ae mn’? © ; 2.2 m = 
ff aed Se? | , or =| Sran{C> 0)= 0 


(1.1) 
Those values Amn! = dan for which equation (1.1) admits solutions that 
are finite at n=+ 1 are the eigenvalues of the differential equation (1. 1). 
The associated eigenfunctions San’? n) are the prolate spheroidal angle 
functions of the first kind, of order m and degree n. Replacement of c 
by-ic in equation (1.1) leads to the oblate spheroidal eigenvalues Amnt"ic) 
and angle functions of the first kind Sa nics n). 

When c vanishes, the differential equation (1.1) becomes the one 
which is satisfied by the associated Legendre functions. It follows that the 
angle functions of the first kind must reduce to the associated Legendre 
functions Pon * of the first kind of integral order and degree, as c 


goes to zero. Therefore 
A (Oh = nin+1), n >m, (1.2) 


When c 1s not zero, equation (1.1) differs from the associated 
Legendre equation by having an essential singularity at infinity. This 
suggests for the angle functions of the first kind, an infinite sum of the 


form oo, 
hae Ae Mmny.y pM 7.\ 
Sane? 7 2. dc) PL dn (1. 3a) 
™) 
Suits m= ; dn(-ic) Pet) (1. 3b) 
r=u, 
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Il NUMERICAL DETERMINATION OF EIGENVALUES 


Eigenvalue Equations 5 Siete 
In his book "'Spheroidal Wave Functions, '' Flammer [1, p. 16 | studies 


the differential equation satisfied by the angle furctions Sant? n}, where 


ONE Tae 


in his notation 


dS__-c, 7) 2 
d 2. mu’? : 22 m ” 
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(1. 1) 
Those values dmnn'£! = on for which equation (1.1) admits solutions that ; a 
are finite at n= +1 are the eigenvalues of the differential equation (1. 1). ona 
The associated eigenfunctions San’? n) are the prolate spheroidal angle 
functions of the first kind, of order m and degree n. Replacement of c 
by-ic in equation (1.1) leads to the oblate spheroidal eigenvalues Amn'ic) 
and angle functions of the first kind S man io n). 

When c vanishes, the differential equation (1.1) becomes the one 
which is satisfied by the associated Legendre functions. It follows that the 
angle functions of the first kind must reduce to the associated Legendre 
functions Pm * of the first kind of integral order and degree, as c 


goes to zero. Therefore 
mn =nin+l), n>m. (1.2) 


When c is not zero, equation (1.1) differs from the associated H 
Legendre equation by having an essential singularity at infinity. This 
suggests for the angle functions of the first kind, an infinite sum of the 
form co 


; 

mn, m 4 , 

Sane? 7 + De, d. (ce) Pn) (1.3a) § 3 
3 


r-0 
Saget m= d"(-ic) Pn) (1. 3b) 
i r= jl 


m m+r 


The prime over the summation sigs: :nd:cates that the summation 1s 
over only even values of r whei.n-m :s €ven, aud over only odd values 7 
of r when n-m is odd. Subsiitution of equation 11. 3a: in equation (1.1) 
with the subsequent use of the associated Legendre diiferential equation 
and of the recursion formulas ior the associated Legedre functions, yields 


the following recursion formula for the coefficients dae 


. Fa . 2 
‘2m+r+2:(2m+r+i1'c” mn, . 


‘2m + 2r +3, (2m + 2r + 5; qn eg 


om Sa sai ea eon a 2 


{ ° ¥ . ae : mn 
+|imer) mere dy Amon! ss (2m +2r- 1i'2m + 2r + 3) © qd. (c) 
+ r(r - 1)c ge 0. tr ene m 0) 44) 
(m+ or-onam+or-1)%r-2 > % Fr a ae 


This recursion formula constitutes a linear homogeneous difference 
equation of the second order. A secona-order difference equation corresponds 
to a differential equation of second order so that there are two non-trivial 
independent solutions. In equation ‘1,1! as r approaches infinity, either 
de / Ba Increases as -4r” / cf or goes to zero as see) (ar). Of 
these two solutions, the former leads to a divergent series; therefore, the 
latter 1s chosen. The expansion (1.3a° ther converges absolutely for all 
finite 7. 

Atter this brief resume of how Flammer aerivea equation (1.4), itis 
rewritten for purposes of numer; a: COmput2uon as follows: 


The difference equation 


'r=n-m 20) (1.5) 


where 
’ , 2 
rr-ljc 


Ais Fe ee at ee ee (1.6) 2 
[2(m + r) - 3 | |2im +r} - 1 | 
, , 2 
t V¢ - - : 
Bo=im+rj‘mer+i)+2mtn mers tj-am =) 2 (1.7) 
[2(m +r}-1][2an+r)+3 | 
[r+2m+2|[r+2m +1 | 9 
ooo (1.8) 
C42 [2 (m +r) + 3] [2 ‘m +r) + 5 | 
has for solutions the convergent infinite continued fraction 
d.42 - A, (1.9) 
~* ey= Bn. * OC.  @ : 
a, dmn'® B42 C44 qs 4 
q.42, 
and the finite continued fraction 
ne ; —— tne): 
r mn rnd” png Ornd 
qd... 
wherein the last term 1s a, OI 44 depending on whether r is even or odd. 
d d 
2 3 


The similarity in prolate and oblate eigenvalues and expansion 


coefficients permits the following section to deal primarily with the prolate 5 
case. These equations can be changed from the prolate to the oblate case 


by the substitution of - ic for c. 


The Parameter n 

In the following discussion Amn) is referred to as \. Consider 
c #0, for when ¢ = 0a special case arises in which equation (1.5) reduces 
to the form )= nin+1). 

It appears that the A's for which the left side of equation (1.5) is 
equal tozero, depend upon the values assigned the narameters m, c and n. 
It will now be shown, however, that the X's depend only on the values of 
parameters m and c and the parity of r = n-m. 

Using the definition: 

r=(n-m)>0, where m>0, 

holding m and c constant, and making use of equation (1.9), equation 
(1.5) may be written 


r-2 _ r r+2 
\- B.-A.) Pee ee ad, (1. 5a) 
A- Bie” a 


qo 
Replacing r by r+2 and making use of equation (1.10), equation 
(1.5) may be written 


x-B .-C Seg EE ROE Ss Sa (1. 5b} 
r+2 r+4dio a : ; 


acer ee ORL CL IOO NT OTE On A 
—_— 


Thus (1.52) and (1. 5b) can be considered equations of the respective forms - 


uv-A C 


——t <0 (1.5¢) 
Py 
uv-A_ C 
r “r+2 =90 (1. 5d) 
u 
where d 
= r-2 1.11 
u=’-B -A_, d. (1. ) 
d 
$e LABEBeee. L (1,12) 


r+2 r+4d r32 = 


From equations (1.6) and (1.8) it can be seen that A. C42 #0. 
Therefore, equations (1. 5c) and (1. 5d) must have identical roots although 
their poles are not identical. 

The conclusion is that given any m and c, equation (1.5) will yield 
a different equation for each new even-valued r, but all of these equations 
will have identical roots. Similarily, given odd-valued r's, equation (1.5) 
will yield a family of equations with identical roots. The shapes of the 
graphs of these equations including the location of poles are altered with 
each new r. 

Figures la - 1d display the effect of varying the value but not the 
parity of r. In each example the value of n was increased by two while 
m and c were held constant. Values of dare plotted along the 
horizontal axis, and corresponding values cf F(\), where F()) is equal 
to the left side of equation (1.5), are plotted along the vertical axis. 
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Figure ld 


Figure 1 — Graphs of Equation 1.5 Given m = 0, ¢ = 6, nas Indicated, 
dX = 0-150 in Incremental Steps of .25 
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F(A) 


Aon (9. 5) 
Figure 2d — Superinposed Figures 2a, 2b, and 2c 


Figure 2 — Graphs of Equation 1.5 Given m = 0, c = 9.5, nas Indicated, 
= 25~—175 in Incremental Steps of .25 
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es Numerical determination of the roots of equation {1.5' 1s usually 

e accomplished through some form of a root-seeking technique (i.e. , 

} a linear fit, parabolic fit, hyperbolic fit). The effectiveness of sucha 
technique is dependent upon the shape of the graph of the funct:on in 
the area of the root. In Figure la, only three roots are discernable 
because the steepness of the curve 1n the neighborhood of each root 
increases rapidly as ) is increased, and the widths of the neighborhoods 
where discontinuities occur decrease as \ increases. Only through 
repeated reduction of the \) increment could the eye detect or the technique 
used locate additional roots. Such reductions soon render a root-seeking 
technique impractical. 

The valuc assigned to n must be carefully chosen when calculating 
the roots of equation (1.5) because while each new n of the same parity 
results in a new equation with identical roots, the graph of each of these 
equations is of a different shape. The relationship between n and the shape 
of the graph of the equation (1.5) for n-m even is such that when n=m the 
shape of the graph of the resulting equation in the region of its first positive 
root is favorably affected, see Figure la; when n=m+2 the region surrovnding 
this new equation's second positive root is favorably affected, see Figure 1b; 
when n=m+4 the region about the equation's third root is affected, see 
Figure lc, etc. Acorresponding relationship exists for (n-m) odd, n=m+1, 
n=m+3,... . 

Studies indicate that for certain combinations of m and c, one value of 
n will fail to influence the shape of the resultant graph in its usual manner. 
In these instances the steepness of the curve within the immediate region 
of the desired root is not decreased enough to permit successful use of a 
root seeking technique. The result is that the technique gives an adjacent 
root as a solution. The elusive root can still be located by recomputing 
equation ‘1 5) replacing n by n-2 or n+2 and solving for the Sa root. 

Figure 2b illustrates a failure when m=0, c=9.5, n=4 (second axis crossing). 
This root can be readily located by recomputing equation 1.5 with n=n-2 
‘Figure 2a', or n=n+2 (Figure 2c:. 
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IV. EXPANSION COEFFICIENTS 


Expansion coefficients a may be obtained in terms of an arbitrary 


coefficient, dy or d,, from equa.ions (1.9) and (1.10) using the correct 
value of )\.__: 
mn 


d d 
ace 1 at 2 (2.1) 
a (Be Gl d, 
r=0,2,4,... 
d d d 
: Tt re} _3 
a, -4| 2, re 7 (2.2) 
r=1, 3,5, : 


The arbitrariness is removed by the adoption of the normalization 
used in Flammer. 
The magnitude of d Fe is fixed by 


r n-m 
oo 2 d 2 
" (-1)° (vr + 2m)! or (-1) (n+m)! 
d SS a (n - m) even, (2.3) 
) oy ary + 2m),d, gi-m =m), +My, 
2 2 2 2 
The magnitude of d, is fixed by 
r-l n-m-1l 
= 2 ar: 2 
d (- 1) (r+2m+1)!o or e (- 1) (n+m +1)! (n= ti)odd 
lr=il, r-1,, ~+2amil d n-mn-m-1,,n+m+41,,’ | 
Zr ( yt ye “1 2 (—>—)! (—5——-)! 
2 2 2 2 
(2.4) 
The expansion coefficients d are thereby completely determined. 
~-4 
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V. THE COMPUTER PROGRAM 


Spheroidal wave runction eigenvalues and expansion coefficients 
were programmed on an IBM 7090 digital computer. The program was 
divided into four separate routines: 

1. Executive Routine 

2. Hyperbolic Subroutine 

3. Ordinate Subroutine 

4. Coefficient Subroutine 
Executive Routine 

This routine has two primary functions; to compute equations (1. 6), 
(1.7), and (1.8) for each set of parameters, and to predict the next eigen- 
value to be computed. 

Increasing c with small constant increments from an initial value 
of zero results in smoothly increasing eigenvalues as shown in Figure 3 
and Figure 4. Successive eigenvalues can thus be predicted with a high 
degree of accuracy. Most root-seeking techniques such as the one used 
in this program need an initial estimate of the root's location, and the 
more accurate the estimate the fewer the calculations required. 

As previously noted, certain combinations of m, c, and n may 
result in the calculation of an adjacent rather than the desired root. 
These undesired solutions are recognized during computation by examination 
of first differences of successive eigenvalues. When such a solution is 
encountered, r is replaced by r + 2, and equation (1.5) is again solved. 
The correct ) is then printed out with an asterisk, which signifies the 
presence of a trouble area. 

Hyperbolic Subroutine 

A hyperbolic root-seeking technique was selected after examination 
of the appearance of the graphs of numerous eigenvalue equations. Para- 
bolic and straight line methods were used with limited success. Eigen- 
values are computed with accuracy to approximately 14 significant digits. 
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Ordinate Subroutine 


Ths routine receives approximations of mn ©) from the hyperbolic 
subroutine, and returns corresponding values of F()). 


Coefficient Subroutine 


and expansion coefficients are computed using equations (1.9) and (1.10). 
Normalization is accomplished and coefficients, correct to approximately 
12 significant digits, are printed in both normalized and unnormalized form. 


The speed of the program is such that the eigenvalue and thirty 
corresponding expansion coefficients for a given m, c, and n can be 


j 
5 
; 
. i L 
Correct eigenvalues are received from the hyperbolic subroutine, | 
computed in less than two seconds. : 
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VI. APPENDIX 
The difficulty in computing eigenvalues for certain combinations of 
m, c, andn was discussed in earlier sections. The following is a list 
of a few such combinations. Correct solutions are given. 


Prolate Eigenvalues: 


dy a7 (18.25) = 4. 8441935628 x 10° 
a, 19 (14-28) = 2.2133274256 x 10° 
ig the) = 2. 9662027320 x 10° 
de yg (20.5) 9 = 5. 8271662786 x 10°2 
Meus: = 7,9153702431 x 10° 
dy ig (21-75) = 4,5767279276 x 10°@ 
dy og (245) = 9. 1471897716 x 10° 
Oblate Eigenvalues: 
Xo 19 (20:75) = - 4. 8870880218 x 10° 
dag (160) = 6. 6192328449 x 10°: 
dy qo (22-25) = - 2. 3406085635 x 10°! 
dy 1g (24:0) = - 6. 7428891070 x 10° 
Ge (1080) . «= 9.3990925921 x 107; 
Meg Albeas) “ve 4. 4862352475 x 10° 
dg 9 (21.25) = 3. 0624599862 x 10°2 
de jg (18.0) = 9, 9902450382 x 10°? 
dy iq (18-8) = 1, 0251012270 x 10° 
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Tables of Numerical Values - Introduction 


The tables contain values for the spheroidal eigenvalues and 
expansion coefficients. 

Calculations made were carried out to sixteen significant digits, 
and table values, which were not rounded, are considered accurate to 
the last digit given. Values are printed in floating arithmetic form; 
thus the number 1. 6980967154 E2 is equal to 1. 6980967154 x 10°. 
Blanks in the tables indicate that the values of the quantities were so 
small they were considered insignificant. 

Where possible, table values have been compared for accuracy 
with previously published values. . 
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